Bright soliton trains of trapped Bose-Einstein condensates 
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We variationally determine the dynamics of bright soliton trains composed of harmonically 
trapped Bose-Einstein condensates with attractive interatomic interactions. In particular, we obtain 
the interaction potential between two solitons. We also discuss the formation of soliton trains due 
to the quantum mechanical phase fluctuations of a one-dimensional condensate. 
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Introduction. — The existence of solitonic solutions is 
a very general feature of nonlinear wave equations. Soli- 
tons, therefore, play an important role in a large num- 
ber of areas in physics such as fluid dynamics, optics 
and high-energy physics. In the case of an atomic Bose- 
Einstein condensate, the macroscopic wave function of 
the condensate obeys the so-called Gross-Pitaevskii equa- 
tion, whose nonlinearity is a result of the interatomic in- 
teractions. Depending on the repulsive or attractive na- 
ture of the interatomic interactions, the Gross-Pitaevskii 
equation allows for either dark, or bright solitons, respec- 
tively. Dark solitons are associated with a notch in the 
amplitude of the macroscopic wave function. Bright soli- 
tons, on the other hand, correspond to localized wave 
packets. The most important feature of both types of 
solitons is that they do not disperse in time due to the 
nonlinearity in the wave equation. The properties of dark 
solitons have been extensively studied theoretically . 
They have also been created experimentally in elongated 
Bose-Einstein condensates p-|ll| . Much less is known of 
bright solitons | fl2|| , which have only very recently been 
created in two experiments with Bose-Einstein conden- 
sates of 7 Li atoms |jl3|Jl4) . 

In the experiment of Strecker et al., soliton trains con- 
sisting of up to 10 solitons have been observed jL4j . More- 
over, it was found that even though the interatomic in- 
teractions are attractive, the neighboring solitons repel 
each other with a force that is dependent on their sepa- 
ration. Here we confirm that the source of this repulsive 
force is a phase difference of tt between two neighboring 
solitons |1| . Physically, this can be understood from the 
fact that the antisymmetric nature of the many-soliton 
wave function prevents the solitons from penetrating each 
other. Although the phase difference explains the repul- 
sive interactions between the solitons, it presents us with 
the problem of understanding why the experiment always 
seems to create soliton trains in which the neighboring 
solitons have a tt phase difference. We also address this 
interesting question here. 

To understand the dynamics of the solitons we use a 
variational approach, in which we describe the wave func- 
tion of the individual solitons as a gaussian and then 
make an appropriate linear superposition of these gaus- 



sians to represent the soliton train. Both the width and 
position of a gaussian can be varied. However, since we 
are dealing with solitons, which by definition do not show 
any dispersion, only the position is allowed to depend on 
time. Using this trial wave function, we derive equations 
of motion for the center-of-mass position and the relative 
distances between the solitons. As expected, the center- 
of-mass motion decouples completely from the relative 
motion. The dynamics of the soliton separations compare 
favorably with the experimentally observed behavior. 

Single soliton dynamics. — We first study the equilib- 
rium properties and dynamics of a single soliton. As in 
the experiments, we consider bright solitons formed from 
a condensate that is confined by a very elongated axially 
symmetric harmonic potential, for which the ratio of the 
radial and the axial trapping frequency obeys uj r /uj z 3> 1. 
We are, therefore, mostly interested in the motion of the 
soliton in the axial direction and can represent the wave 
function of a single soliton by a trial function of the form 
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where £(t) represents the soliton's center of mass, and m 
is the atomic mass. For the amplitude of the soliton, we 
take the gaussian 



A so i(x) = 



N 



exp (-^r^)' (2) 



where q z and q r are two variational parameters that de- 
termine the width of the soliton in the axial and radial 
directions, respectively. The total number of atoms in 
the condensate is denoted by N. Note that this gaussian 
ansatz is particularly well-suited when the axial width of 
the soliton is comparable to the axial harmonic oscilla- 
tor length l z — ^h/muj z . In the case that q z <C £ z , it 
would be better to replace the exponential exp(— z 2 /2q 2 ) 
by y / 7r/[2cosh(z/g 2 )] jn|[L3|. However, we will see later 
on that the gaussian ansatz gives physically reasonable 
results in this case also. 

The equilibrium widths of the soliton can be calculated 
by minimizing the Gross-Pitaevskii energy functional 
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with respect to q z and q r Jig ]. The external poten- 
tial obeys U cxt (x) = m(u 2 z 2 + u> 2 r 2 )/2 and the inter- 
atomic interaction strength T 2B = Anah /m is propor- 
tional to the negative s-wave scattering length a. As- 
suming that q r has been found in the above manner, the 
physics becomes effectively one dimensional. The one- 
dimensional energy for ip(z,t) = J dr2irrip(x,t) is again 
given by Eq. (||) but now we must use T 2B = Airnh 2 /m 
with k — a/2irq 2 . It should be kept in mind, however, 
that the one-dimensional theory can only be applied for 
a soliton containing less than iV max = 0{t r /\a\) atoms, 
where i r = W%J muj r is the radial harmonic oscillator 
length. Above this number of atoms the soliton will col- 
lapse 

Using this effective one-dimensional picture, the equa- 
tion of motion for is finally derived by substituting 
the gaussian ansatz ip(z,t) into the action S[ip*,ip] — 
J dt(J dziip* dip / dt - E [ip*, ip]). The resulting equation 
d 2 ((t)/dt 2 = —uj 2 C,(t) shows that the center of mass 
of the soliton oscillates sinusoidially with the trap fre- 
quency. This is, of course, an exact result that follows 
from the Kohn theorem p(J . 

Two soliton dynamics. — Our trial wave functions for 
two solitons with a phase difference of or 7r are, respec- 
tively, the symmetric and antisymmetric combinations 
of the two gaussian wave functions of a single soliton. In 
detail we use 
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A sol (z - Zl {t)) exp ( i-^±l(z -at)) 
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Here Z\ it) and z-i (t) denote the positions of the two soli- 
tons and = [z\(t) + z 2 (2)]/2 is their center of mass. 
For simplicity we consider here only the case of a phase 
difference of or n between the solitons. In addition, we 
take the number of atoms in each soliton, and therefore 
also their widths, equal to each other. The generalization 
to an arbitrary phase and an arbitrary number of atoms 
is straightforward but somewhat tedious. Fortunately, 
the above assumptions turn out to be reasonable when 
we compare our results with experiments, and it allows 
us to bring out the dynamics of the relative separation 
most clearly Finally, for the rest of this section we scale 
length to 



- y^h/muz, time to l/u> z , and energy to 
The normalization constant N± (t) is then given by 



N±(t) = 2 Izbexp 



v 2 (t) 
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where rj(t) — z\(t) — z^it) is the distance between the two 
solitons. Note that N is still the total number of atoms 



in the condensate, so a single soliton contains only N/2 
atoms. 

Using this trial wave function, the lagrangian per atom 
takes the form 



L[C,y] _ l 
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where the effective potential V[rj\ is given by 



V[rj\ = 

Tj^[(l + q^W + (qi-W 2 ] exp 
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and f) = dr)/dt. It is clear from this lagrangian that the 
center-of-mass motion and the relative motion are com- 
pletely decoupled. Furthermore, the equation of motion 
for the center of mass is simply d 2 ((t)/dt 2 = — which 
corresponds again to the Kohn mode that oscillates at the 
frequency of the trap. On the other hand, the equation of 
motion for the relative motion is rather complex due to 
the coupling between the two solitons that is determined 
by the velocity-dependent potential V[r]]. 

In Fig. |l|, we compare the exact numerical solution of 
the relative motion with experimental data. The param- 
eters of the experiment and the theory are the same as 
those of Ref. Q] , and we take the anti-symmetric com- 
bination of gaussians. Very good agreement is obtained 
with both the experimentally observed amplitude of oscil- 
lation of the relative coordinate and its frequency, which 
is approximately 2uj z . In particular, we see that the two 
solitons never cross each other. This result, and the non- 
sinusoidal nature of rj(t), are both due to the hard-core 
nature of the interaction when the solitons start to over- 
lap. 

To show this, we neglect the velocity dependence of 
V[rj\, assume that q z <C 1 and take the limit of large dis- 
tances between the two solitons, i.e., rj ^> q z . With these 
approximations we obtain 



r . 2V2ttNk / 4V27rA^K 
V[rj\ = T 



H exp 
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The second term in the right-hand side represents the 
interaction energy V(r]) between the two solitons. For 



2 



the antisymmetric choice of the two-soliton wave func- 
tion this potential energy is positive, which implies that 
the two solitons will indeed repel each other. 




FIG. 1. The relative motion of two solitons. There are 5000 
atoms in each soliton and a = — 3ao, with ao the Bohr radius. 
The solid curve is the solution of the exact equation of mo- 
tion for n(t). The dashed curve is the solution of Eq. (|^) with 
V(n) = (n 2 /2qf) exp( — n 2 /<±q 2 ). The points are experimental 
data obtained using the same apparatus described in Ref. [14]. 
In the experiment more than two solitons were created. We 
focused on two large and adjacent ones. 

In this approximation, the equation of motion of r](t) 
takes the form of Newton's equation 

<?g) = (4) _ mm (9) 

dt z drj 

Thus, the problem reduces to that of a particle moving in 
the potential rj 2 /2 + V(rf), which has two minima located 
symmetrically around 77 = 0. Depending on the initial 
conditions, the two solitons can be trapped in these min- 
ima. In other words, the separation between them oscil- 
lates and the solitons never cross each other, as verified 
by the exact solution. To investigate the importance of 
the interatomic interactions, we show in Fig. [l] with the 
dashed curve, the solution to Eq. when the potential 
contains only the contribution from the kinetic energy of 
the atoms. A solution of Eq. (0) with the full potential 
V(rj) from Eq. (^) is almost identical to the exact solu- 
tion. The smaller amplitude of oscillation in the relative 
motion is, therefore, a clear signature of the reduction 
of the repulsive force between the solitons due to the 
interatomic interactions. We have also analytically cal- 
culated the interaction energy for two solitons with an 
amplitude proportional to y/w/[2cosh(z/q z )] and solved 
the corresponding equation of motion. It leads to only 
minor corrections that are due to the fact that the po- 
tential V(ri) is now somewhat more repulsive and only 
falls off as (4?j/g|) exp(— rj/q z ) for large separations. 

Soliton train formation - - In the experiment of 
Strecker et al. Q, soliton trains are formed by making 



use of a so-called Feshbach resonance [|2l| . In summary, 
the procedure consists of first making a stable condensate 
with a relatively large positive scattering length, and then 
switching to a small negative value. An important clue 
for understanding the mechanism for the creation of the 
soliton train is the observation that all the solitons repel 
each other, and thus have a phase difference equal to n. 
This suggests that by changing the sign of the scattering 
length, the phase modes of the condensate become unsta- 
ble. Roughly speaking, the most important wavelength 
of the unstable phase modes must then be comparable to 
the separation between solitons in the train. 

To make this physical idea quantitative, we consider a 
homogeneous condensate with a density n — N/L, where 
L is the (Thomas-Fermi) length of the condensate. De- 
noting the phase of the condensate wave function ip(z, t) 
by x(z,t), it is easy to show from the action S[ip*, ip] 
that Xk(t) = / dz e~ tkz x(z,t) corresponds to a (com- 
plex) harmonic oscillator for every wave vector k > 0. 
Physically this implies that the quantum mechanics of 
every phase mode is analogous to the quantum mechan- 
ics of a fictitious particle in a two-dimensional harmonic 
oscillator potential. The mass of this fictitious parti- 
cle is mk = Ae]~n/uj'j,L and the frequency of the har- 
monic oscillator potential is uj k — y € \ + 2nT 2B e k /%, 
with e k = H 2 k 2 /2m. Initially the scattering length is 
positive and all these frequencies are positive. Assum- 
ing that the condensate has come to equilibrium, the 
quantum fluctuations of Xk(t) are then fully determined 
at zero temperature by the ground-state wave function 
(l/ir£l)exp(-\xk\ 2 /2£ 2 k ), with 4 = yfh/m k uj k . In par- 
ticular, we have that (x{z, t)) = 0. 

However, if at t — we suddenly change the scattering 
length to a negative value, the modes with k < k max — 
y/l6w\n\n become unstable because the spring constant 
m k 0J k of these harmonic oscillators becomes negative and 
equal to — m k fl k . Solving the quantum mechanics of a 
particle in an inverted harmonic oscillator potential, we 
can then show that for t S> l/f2fc 

<l«< f )l 2 > = 4^l( 1 + C e2 °"- <10 ' 

As in the case of a spontaneously broken symmetry, the 
latter result can physically be understood by saying that 
quantum mechanical fluctuations imprint the condensate 
wave function with the phase 

(11) 

Having obtained this result, we are now able to sim- 
ulate the experiments. We first determine the ground 
state wave function of a condensate with a positive scat- 
tering length for the experimental parameters of interest. 
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We then imprint this wave function with the phase given 
in Eq. ([ll]) and let the resulting wave function evolve 
under the appropriate Gross-Pitaevskii equation with a 
negative scattering length. We note here that, due to 
numerical limitations, we have used an initial number 
of 10 4 atoms which is an order of magnitude less than 
that in the experiment. In this case, we find that af- 
ter an evolution time of t ~ 1.8/uj z , seven solitons are 
formed, together with phonon excitations. The solitons 
are clearly visible as plateaus in the phase, after subtrac- 
tion of the parabolic behaviour of the phase due to the 
overall shrinking of the condensate. The phase difference 
between two adjacent solitons is indeed approximately 
equal to 7T. Some deviations from ir are expected because 
the solitons are moving relative to each other. These re- 
sults are shown in Fig. ||, where we plot both the density 
and the phase of the condensate wave function. Note 
that our calculation includes no relaxation processes due 
to the presence of a thermal cloud |p2[ . Including these 
stabilizes the soliton train by damping the phonon exci- 
tations and the soliton motion. 
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FIG. 2. Soliton train formation. The solid curve is the 
density and the dashed curve is the phase of the condensate. 
Initially, we start with an equilibrium condensate profile for 
10 4 atoms with a scattering length of 200ao- The scattering 
length is then instantaneously changed to a = — 3ao. The 
trap parameters are obtained from Strecker et al. [14]. 

Conclusions — In this Letter, we have considered the 
formation and dynamics of bright soliton trains in a 
trapped Bose-Einstein condensate with attractive inter- 
actions. In particular, we have shown that two solitons 
with a phase difference of ir repel each other. The result- 
ing dynamics compares favorably with experiment. In 
the experiment, always more than two solitons were cre- 
ated, while here we have focused on only two. In princi- 
ple, a generalization to more than two solitons is straight- 
forward. In this case, we just sum in the lagrangian of 
Eq. (^) over all the positions Zj(t) of the bright solitons. 

We have also shown that soliton trains can be dynam- 
ically generated by quantum mechanical phase fluctua- 
tions of the condensate, due to the modulational insta- 



bility [g3[ that exists in a Bose-Einstein condensate with 
attractive interactions. This mechanism gives an intu- 
itive explanation why adjacent solitons have a n phase 
difference. It also predicts that the number of solitons 
that is formed after a rapid, i.e., fast compared to l/w z , 
sign change in the scattering length, is equal to the ra- 
tio of the initial condensate size and the wavelength of 
the most unstable phase mode. However, this prediction 
does not take dissipation into account. Because solitons 
are not topological objects, dissipation can lead to a con- 
siderably different number of solitons as it drives the gas 
to equilibrium. In the experiment, the number of solitons 
in the train is about half of the dissipationless prediction 
. This observation suggests that dissipation processes 
are indeed present. A detailed comparison with exper- 
iment thus requires knowledge of the properties of the 
thermal cloud, which are not well known experimentally. 
It was also observed that the number of solitons in the 
train is, over a large range, independent of the expo- 
nential time constant with which the scattering length is 
changed. This is as expected as long as dissipation plays 
an important role and the sign change in the scattering 
length remains nonadiabatic. 
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